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(*C) Abstract 

We study the 3D forced-dissipated Gross-Pitaevskii equation. We force at 
relatively low wave numbers, expecting to observe a direct energy cascade 
^ and a consequent power-law spectrum of the form k~ a . Our numerical results 

show that the exponent a strongly depends on how the inverse particle cas- 
cade is attenuated at fc's lower than the forcing wave number. If the inverse 
1—1 cascade is arrested by a friction at low fc's, we observe an exponent which 

i is in good agreement with the weak wave turbulence prediction k~ x . For a 

hypo-viscosity, a k~ 2 spectrum is observed which we explain using a critical 
balance argument. In simulations without any low- A: dissipation, a conden- 
sate at k = is growing and the system goes through a strongly-turbulent 
transition from a four-wave to a three-wave weak turbulence acoustic regime 
with /c -3 / 2 Zakharov-Sagdeev spectrum. In this regime, we also observe a 
spectrum for the incompressible kinetic energy which formally resembles the 
Kolmogorov fc _5//3 , but whose correct explanation should be in terms of the 
Kelvin wave turbulence. The probability density functions for the velocities 
^ and the densities are also discussed. 

C3 Keywords: 



1. INTRODUCTION 

After the seminal papers by A. Kolmogorov in 1941, it is well established 
that, apart from small corrections due to intermittency p], the energy spec- 
trum, E(k), of the velocity fluctuations for high Reynolds number hydrody- 
namic turbulence shows a power law of the form E(k) = C P 2 / 3 k~ 5 ^ 3 , where 
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C is the dimensionless Kolmogorov constant and P is the flux of energy in 
the wave number space. This is a very strong result that has been confirmed 
experimentally and numerically by the direct numerical integration of the 
Navier-Stokes equation. It can be obtained via dimensional considerations 
or as a solution of phenomenological turbulence closures pQ . However, so far, 
this result has not been obtained analytically from the Navier-Stokes equa- 
tion. Many decades after the work by Kolmogorov, it has been discovered 
by Zakharov in 1965 that systems of weakly nonlinear, dispersive, random 
waves behave qualitatively in a similar way as hydrodynamical turbulence 
[2]. Namely, the nonlinear interaction of waves can produce other waves 
with different wavelengths and so on, generating a cascade process leading 
to power-law wave spectra similar to the Kolmogorov spectrum. Because of 
such an analogy, they are called Kolmogorov- Zakharov (KZ) spectra, and 
the entire field is known as Weak Wave Turbulence (WWT). Description of 
WWT turns out to be more accessible than of the hydrodynamic turbulence 
because nonlinearity in the dynamical equations, although still crucial, is 
small. In the WWT framework, a systematic approach based on averaging 
the dynamical equations leads to a Boltzmann-like equation known as the 
wave kinetic equation which describes the evolution of the spectrum of the 
turbulent wave field [3]. One remarkable property of WWT is that, in con- 
trast to hydrodynamic turbulence, the KZ spectra have been found as exact 
stationary solutions of the wave kinetic equation [2JH]. Unlike the thermody- 
namic solutions, for which the integrand in the collision integral is identically 
zero, the KZ solutions correspond to non-trivial states for which a source, a 
sink and a window of transparency (inertial range) are required. Since this 
discovery, weak wave turbulence has found applications for a vast variety of 
physical systems ranging from quantum to astrophysical scales, see books 
[3j E] and references therein. 

In this paper, we consider nonlinear dispersive waves described by Gross- 
Pitaevskii equation (GPE), 

<^ + vV + ffMV = o. (i) 

This partial differential equation has attracted the attention of many re- 
searchers: it describes propagation of optical pulses in nonlinear media (6] 
and weakly interacting boson gases at very low temperatures called Bose- 
Einstein condensates (BEC) [7]. In the present paper we will be concerned 
with the latter case and we will focus on the three-dimensional (3D) systems. 
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Complex wave function tp is called order parameter and a = ±1, depending 
on the physics of the problem: the defocusing case, a — — 1, represents re- 
pulsive bosons, while the focusing, a — 1, considers an attractive interaction. 
GPE has recently attracted the attention of many fluid dynamicists because 
BEC is a good example of superffuid, i.e. a fluid with zero viscosity. Indeed, 
using the Madelung transformation, the GPE can be mapped onto an Euler 
equation which differs from the classical one only by an extra term named 
quantum pressure. Therefore, numerical computations of the 3D defocusing 
GPE have become a tool for investigating superfluids and quantum fluids 
dynamics. Phenomena such as the vortex reconnections [S], formation of a 
condensate [9], and formation of power-law spectra [TUJ [TTJ [T2J [32] have been 
observed. 

The purpose of this paper is to revisit and investigate the turbulence char- 
acterises in the defocusing 3D GPE with particular attention to the forced 
and dissipated case. An interesting issue to be addressed is verification of the 
WWT theory which offers a solid theoretical tool for predicting statistical 
quantities in systems of dispersive, weakly interacting waves. As it will be 
shown in the next section, when the nonlinearity becomes large, the predic- 
tions of the WWT theory fail and the concept of critical balance (CB) has 
to be introduced in order to explain the observed spectra. 

This paper is organized as follows. In section [2] we present the theoretical 
background on the GPE model revisiting some mathematical aspects, includ- 
ing quantum vortices and the general properties of quantum turbulence. In 
section [3] we introduce the forced-dissipated GPE and present predictions for 
steady turbulent states: in particular we discuss the WWT for the four-wave 
and the three- wave regimes and the CB conjecture. Section [4] is dedicated to 
presenting numerical results in three different regimes: free condensation at 
large scales (RUN 1), dissipation by friction at low wave numbers (RUN 2), 
and dissipation by hypo- viscosity at low wave numbers (RUN 3). Finally, 
section |5] contains the conclusions. 



2. Theoretical background 

In this paper, we will consider the defocusing GPE model, 

ld ^r~ + v2 ^ x '*) - i^)i 2 ^(x,t) = o, (2) 

where the nonlinearity comes from the self-interactions proportional to the 
gas density p(x, t) = ^(x, t)| 2 ; this term is a consequence of considering local 
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interactions between bosons. The system is conservative and its Hamiltonian 
is 

H = J ^\V^,t)\ 2 + ^(x,t)\^dx = H lm (t) + H nl (t). (3) 

In the latter relation the total energy H has been split into a part responsible 
for the linear dynamics, H[ in {t) = f |Vt/>(x, £)| 2 <ik, and the one describing the 
nonlinear interactions, H n i(t) = 1/2 j |?/>(x, £)| 4 <ik. In the following, we will 
use the energy densities defined as £u n = Hu n /V and £ n i = H n i/V where V 
is the total volume. The system conserves a second quantity, the mass M, 
defined as 

M = J |^(x,t)| 2 rfx = J p(x,t)dx. (4) 

An important quantity that characterizes the system is the healing length 
£: it physically estimates the distance over which the field ^(x, t) recovers its 
bulk value when subject to a localized perturbation. This definition refers to 
the case of a single perturbation in a uniform field but can be extended to 
even highly perturbed statistical systems as 

e = -ijn («) 
v(p) 



where (•) denotes the spatial averaging, i.e. (p) = M/V. The healing length 
£ measures on average the scale at which the nonlinear term becomes compa- 
rable with the linear one. Dual to the scale £ is wave number = 2ir^ (p). 

The GPE has been widely studied in the fluid dynamics framework. In- 
deed, the Madelung transformation, t) = p(x, t)e %e<yyi,t \ maps the GPE 
for the complex field into the system of two equations, 

| + V-(pv) = 



\ dt dxk J dxj dxk 

for the real density field p(x, t) and a real velocity field v(x,i) = 2V0(x,i). 
The first equation represents a continuity equation for a compressible fluid 
and the second equation is a momentum conservation law. The terms in 
the r.h.s of the latter equation can be thought as a pressure p = p 2 and a 
"quantum stress" tensor, Sjfc = P aJ.g^ • The quantum stress term becomes 
important at scales of the order of £. The system (pi) describes an inviscid 
and irrotational fluid flow. 
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Figure 2: p(x, t) near a 2D vor- Figure 3: Plot of the low density region p(x, t) below 
tex. a certain threshold in a freely decay simulation. 



2.1. Quantum vortices 

Even if the fluid is irrotational, particular types of vortex solutions exist. 
This is true if the region occupied by the irrotational flow is not simply 
connected, e.g., if there are phase defects in the field at locations of the 
zero density. Moreover, differently from classical fluids, such vortices carry 
quantized circulation and therefore are called quantum vortices. To better 
understand their structure, we first consider a two-dimensional (2D) system. 
A necessary condition to assure the continuity of the complex field ip is that 
the phase changes by A8 = 2irn, where n G A/", around a vortex. As the 
velocity field v is proportional to the gradient of the phase, it is easy to see 
that the circulation around a vortex is quantized, 

C= <j)v-dl = 2 <j)V6-dl = 2A6. (7) 

In Fig. [I] and Fig. [2] we show respectively the phase field 6*(x, t) and the 
density field p(x, t) in the neighborhood of a 2D vortex embedded in a uniform 
density field. The vortex core size is of the order of £ because, by definition, 
the healing length measures the size of a generic order-one fluctuation on 
the uniform condensate solution. By measuring circulation it is possible 
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to distinguish between clockwise and anti-clockwise vortices. Note another 
important difference from the classical fluids arising from the presence of the 
quantum stress: vortices with the opposite sign can approach each other and 
annihilate. 

In 3D, the vortices are more complicated objects consisting of continu- 
ous lines or loops with different topologies. These structures can oscillate 
producing Kelvin waves [TS] , be transported by the fluid or induce a fluid 
motion, and reconnect [HQS]. The vortex energy is proportional, in the lead- 
ing order, to its length and it can be transferred to the fluid via sound waves 
|16j . Examples of all these vortex motions is shown in Fig. |3j here we plot 
a snapshot of low density regions, associated with the vortex cores, obtained 
in a numerical computation of the GPE in a periodic cubic box with twelve 
vortices as initial condition |17j . 

2.2. Quantum Turbulence 

Vortex-vortex dynamics (collective dynamics of vortex bundles, recon- 
nections of vortex lines), vortex-sound (radiation and scattering of sound 
by vortices) and sound-sound interactions (acoustic wave turbulence) and 
the dynamics of the vortex itself (Kelvin waves, vortex rings) represent the 
ingredients of quantum turbulence (QT). When forced at the scales larger 
than the mean inter-vortex separation £ 3> £, QT cascade starts as the 
classical Richardson cascade characteristic of Navier-Stokes turbulence un- 
til it reaches the scale £. At scales < £ the discreteness of the quantized 
vortex field becomes essential, and the further cascade to lower scales is 
very different than the one in the classical fluids with continuous vorticity 
[T5j [T5j I2U] . After some reconnections and other crossover processes near the 
transitional scale £ [ZI1E2], ^ ne ener gy cascade is be carried to the scales < £ 
by Kelvin wave cascade until, at a very low scale, it is radiated into sound 

[ISl HH1 EQl E31 EH ESI Eg EZ!- 

In the recent years, GPE has become a popular model for studying QT 
by using numerical simulations [TpJ ED DEI EE1 EH] • The pioneering paper 
of Nore et al. [TD] showed similarities between the 3D GPE and the classical 
Navier-Stokes turbulence, including observations of the famous Kolmogorov 
5/3-spectrum. This spectrum was observed for scales > i (with £ > £), which 
agrees with the view that at this scales the discreteness of the vortex lines 
is inessential, and the turbulence picture is basically classical Kolmogorov. 
However, in several follow-up works [TTl IT2| EE] the 5/3-spectrum was re- 
ported to extend to the smaller scales, between £ and £, where the classical 
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Kolmogorov picture is not expected to be valid theoretically. Later in our 
paper we will comment on this apparent paradox. 

Numerical quantities usually measured in QT can be defined in the fluid 
dynamic framework. Namely, the total energy is divided into kinetic, quan- 
tum and internal energy of the system as follows [TU] , 

H = E kin {t) + E qu (t) + E int (t), (8) 

where Eki n {t) = \ J (-^/pv) 2 <ix is the kinetic energy , E qu (t) = | J (2V ^fp) 2 (ix 
is the quantum energy and E int {t) = \ J p 2 d~x = H n i(t) the internal energy. 
Moreover, the quantity in the integrand of the kinetic energy is usually di- 
vided, using the Helmholtz's theorem, into a solenoidal (incompressible) and 
an irrotational (compressible) vector fields, 

v/pv = (y/pv)* + (^pv) c . (9) 

The incompressible part, satisfying V ■ (^/pv) ! = 0, is associated to the 
vortex dynamics; the compressible part, which satisfies V x (y^v) = 0, is 
related to the sound waves. Note that the mentioned above results of papers 
[Till [TT| [T2l [28] for the fc~ 5//3 power-law behavior refer to the spectrum of the 
incompressible kinetic energy E l kin . 



3. Forced-dissipated GPE 

The aim of this research is to study the steady turbulent states in the 
GPE model, with particular interest in observing cascades from large scales 
to small ones. Therefore, in the spirit of classical turbulence, we build a 
sustained system by including a forcing and a dissipation terms in the GPE: 

ld ^&r~ + v2 ^ (x ' t] ~ '^ (x ' t] |2 ^ (x ' t] = T + p ' (10) 

where the forcing T injects mass and energy, while the dissipation T> removes 
them. As a consequence, the mass and the energy are no longer constants 
of motion. However, if a steady turbulent state is reached, the quantities 
H(t) and M(t) are zero. In such a steady state, the mass and the energy get 
transferred from the forcing to the dissipation scales, being approximately 



conserved while cascading though the inertial ranges (see Appendix A) 
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3.1. Weak wave turbulence predictions 

The WWT furnishes some prediction for wave spectra in steady states. 
In order to introduce WWT, let us put the system into a 3D periodic box 
and write the equation (neglecting T and V) in Fourier space: 

% d ^t) _ u{ y x) ^ t) = ^(k 2 ,t)^(k 3 ,t)^(k 4 ,t) ( 5(k 1 +k 2 -k3-k 4 

k 2 ,k3,k 4 eZ 3 

Here ipi = ip(\ti,t) is the Fourier transform of ^(x, t) and u(k) = k 2 is the 
dispersion relation of the system (hereafter k = |k|). The WWT theory, 
taking the limits of infinite-box limit, small nonlinearity, and assuming space 
homogeneity and random phases and amplitudes (RPA) of the initial wave 
amplitudes, provides a statistical closure. The closure predicts the behav- 
ior of statistical quantities such as correlators {iftifa ■■ ■i>ni>n+ii , n+2- • -^n+m) > 
where the average is performed over the random initial data. In this frame- 
work the simplest non-trivial correlator is the second order one, 

(ftfa) = (|&||&|e'<*-»'>) = n^i - kj), (12) 

where the quantity m = n(ki,t) is called wave-action spectrum. The RPA 
assumption allows one to close the system by the Wick decomposition, a 
mechanism which splits the higher order correlators as sums of products of 
second order correlators [30l EU [32] . This procedure, when applied to the 
GPE model, leads to the following four-wave kinetic equation for the evolution 
of the wave-action spectrum [33] . 

dm f ( 1 1 1 1 

47r / m^n^m 



dt J \m n-2 n 3 n 4 ^ 

x5(ki + k 2 - k 3 - k 4 )5(o;i + uo 2 - - ^4)dk 234 , (13) 

where dk 234 = dk 2 <ik 3 <ik 4 . 

This is an integro-differential equation which models "wave collisions" , in 
analogy with Boltzmann kinetic equation for the particle collisions. Physi- 
cally, it says that for times much longer that the fast wave periods Tj = 27r/wj, 
the wave amplitudes ipi are effectively coupled only if they satisfy the follow- 
ing resonant conditions, 

ki + k 2 = k 3 + k 4 

(14) 

Ui + u 2 = OJ3 + OJ4. 
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Equation ( 13 ) has the following invariants, M = J n(k, t)dk, P = J kn(k, t)dk, 
and E = Ja7(k)n(k, t)dk that are respectively the total number of particles 
(or the mass), the momentum and the energy. Moreover, the dynamics de- 
scribed by the kinetic equation is irreversible in time and an entropy measure 
can be defined. Three trivial functions that cancel the integrand in (13) are 
n W = ^ n P) _ ( B . k ^-i and n (3) _ Ck -2 with A B and q constants _ 

By combining these solutions, the general thermodynamic solution takes the 
form: 

^ = ^+(u-k) + ife2' (15) 

where /i is a chemical potential constant, u is a constant macroscopic velocity 
and T is a temperature of the system. In an isotropic field, the macroscopic 
velocity is zero and the relation (15) assumes the form known in literature 
as Rayleigh- Jeans (RJ) thermodynamic equilibrium distribution. As the dis- 
tribution is not convergent for large fe's, the temperature and the chemical 
potential can be evaluated based on the known mass and energy of the system 
only by introducing an ultraviolet cutoff, as proposed in [34J. 

In presence of an external forcing and damping, besides the RJ distribu- 
tion, two other steady solutions may exist, namely the KZ spectra discussed 
in the introduction. These solutions have the form of a power-law n& = ck~ a , 
where c is a dimensional constant. They can be obtained analytically by 
applying a change of integration variables known as the Zakharov transfor- 
mation [3], or by a dimensional analysis [35J. The value of the exponent a 
depends on the particular wave model. KZ solutions correspond to constant 
fluxes of positive conserved quantities in the scale space (turbulent cascades). 

The GPE model has two positive invariants, the mass M and the energy 



E. It can be shown (see Appendix A) that energy has a direct cascade, 



i.e. from large to small scales, while the mass cascade is inverse, from small 
to large scales. The corresponding KZ exponents are respectively = 3 
and olm = 7/3. Hereafter in this work, we present results in terms of the 
one-dimensional spectrum nn)(k), obtained from n(k) by integrating out the 
angular variables, i.e. ni£>(k) = 47r/c 2 rifc. In this notation the KZ solutions 
are 

n[ E I l{k)=AKc^k- 1 

n[ M D \k)=A^k- l '\ 1 ' 

These solutions represent respectively constant fluxes of M and E and, there- 
fore, they are sustained by an external forcing and a dissipation. Even tough 
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the KZ spectra are found for infinite inertial ranges, they are expected in 
finite systems provided the scales of J 7 and D are widely separated in Fourier 
space and provided the interactions of the wave modes are local. The locality 



of the KZ solutions ( 16 ) is checked in Appendix B the energy cascade turns 



out to be marginally nonlocal and its locality is restored by a logarithmic 
correction, while the inverse particle cascade is local. 

3.2. Transition to three-wave interactions 



Equation (13) describes a four- wave interaction process which is respon- 
sible for the direct and the inverse cascades. When the inverse cascade is 
not damped at low fc's, it leads to accumulation of particles at these scales 
which can alter the four- wave dynamics. Respectively, the zero-mode n^^it) 
(related to the uniform part of the field ip in physical space) will grow, 
which can be interpreted as a Bose-Einstein condensation process. When 



the condensate fraction becomes large, the kinetic equation (13) ceases to be 
valid. Suppose that a large fraction of the wave-action is present at the zero- 
mode, thereby ^(x, t) = c(t) + e0(x, t), where represents small fluctuations 
(e < 1). By substituting this ansatz into GPE ({!]), we find, at the order e°, 
the evolution equation for the condensate fraction, 

*t-|c| 2 c = 0. (17) 



dt 



Its solution is c(t) = coe~ ,c o*, where cq is a real positive constant. Thus, the 
condensate amplitude rotates in the complex plane with an angular velocity 
proportional to its square modulus. In the next order in e, we obtain a linear 
equation for the fluctuations on the condensate background: 



«^(x,*) + V 2 0(x,t) - 2c 2 1 0(x,t) + c e- 2iC V(x,t) = 0. (18) 

Diagonalizing the linear dynamics in Fourier space, one can show that the 
linear wave modes oscillate at the Bogoliubov frequency 



w(k) = ±k Jk 2 + 2c 2 . (19) 

In the limit of small k's or strong condensate fraction Co, when k 2 -C C 2 ,, the 
fluctuations are acoustic waves with the speed of sound oj/k ~ a/2cq. 
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In the next order in e, when weak nonlinearity is taken into account for 
the fluctuations, it is possible to use the WWT theory and derive a kinetic 
equation describing three- wave interactions of the Bogoliubov sound [33j EH] : 

dn^ r 

-Qf= (^231 - #123 - R312)dk 12 , (20) 

where 

#123 = 27r|Vi23| 2 <$(ki - k 2 - k 3 )8(ux - cj 2 - u 3 )(n 2 n 3 - n x n 2 - n x n 3 ) (21) 
and the analytical form of the scattering matrix V123 is given in [36]. Equation 



(20) describes three- wave interactions 1^2 which conserve only the energy 
and not the mass. Thus, only the energy cascade KZ solution is relevant 
to this regime. For the large-scale (strong condensate) limit k 2 <C Cq, the 
direct cascade KZ spectrum takes the form E\r>{k) ~ &;~ 3 / 2 predicted by 
Zakahrov and Sagdeev for the 3D acoustic WWT [37J. Because most of 
the wave-action in this case is in the condensate, the ID energy spectrum is 
Eu)(k) ~ c^n^](k) [3S]- Therefore, the wave-action spectrum for the energy 
cascade in the acoustic regime is 



n 



(E) 
ID 



(k) ~ AT 3/2 . (22) 



3.3. Critical balance conjecture 

In some physical situations, wave turbulence fails to be weak, and the 
wave spectrum saturates at a critical shape such that the linear term is 
of the same size as the nonlinear term for each mode k. Such a critical 
balance appears to be typical for a wide range of physical systems, ranging 
from Magneto-Hydrodynamic turbulence |38j, to the rotating and stratified 
geophysical systems [39] . The most famous example is the Phillips spectrum 
of the gravity water waves [IH1 SI] , in which case the saturation at the critical 
value occurs due to wave breaking. 

In the GPE model, similar situation may occur when an equivalent of 
wave breaking process is active in the system. Namely, we will see that when 
the low-k range is over-dissipated by strong hypo- viscosity (RUN 3), the 
inverse particle cascade tends to accumulate at low k's (infrared bottleneck) 
until a critical balance is reached and the spectrum is saturated. Indeed, 
for the inverse cascade to exist the wave turbulence must be weak, because 



only then the Fj0rtoft argument works (see Appendix C). However, when the 
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linear and the nonlinear terms, locally in Fourier space, are of the same order, 
the inverse cascade stops and so does the infrared bottleneck accumulation. 
This is precisely the mechanism of reaching the critical balance condition in 
this case. 

Now we present an estimate for the critical balance spectrum in the GPE 
model. Equating the linear and the nonlinear terms in Fourier space gives 

=► k-* ~ |^| 2 . {26} 

Note that we have replaced each dk integration by k 3 thereby assuming that 
only the wave amplitudes with similar fc's are correlated in Fourier space. 
Thus, for the ID wave-action spectrum in the critical balance regime, we 
have 

nf D B \k) ~ AT 3 . (24) 



4. The numerical experiments 

As we want to understand the basic properties of the GPE turbulence, 
we choose to deal with the simplest configuration: triple periodic boundary 
condition and uniform mesh grid. With this choice we can use the discrete 
Fourier transforms which are numerically fast (421 143] . In our numerics the 
complex wave field ip is a double precision variable defined in space over 
256 3 points. The simulation box has side L = 256 and so the Fourier space 
width is Lfc = 2tt with resolution AA; = 27r/256. Without considering for the 
moment the forcing and the dissipation terms, the GPE model ([!]) can be 
written in the physical space as a sum of a linear operator £ and a nonlinear 
one Q 

i^ = C^ + G^, (25) 

where C = —V 2 and Q = \ip\ 2 . We then use a split step method to solve 
separately the contributions of the two operators in time. This choice is 
very useful because the linear operator has the exact solution in Fourier 
space, 4>(k,t + At) = i/>(k, t)e~*' k ' At , while the nonlinear one, due to the 
conservation of the mass in the system, has the analytic solution ip(ya,t + 
At) = ^(x, t)e~*'^ ( ' x '*- ) ' At in the physical space. At each time step At, we 
first evaluate the linear part in Fourier space and then use this temporary 
solution to solve the nonlinear part in the physical space. With this choice, 
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cases 


forcing 


dissipation at low fc's 


dissipation at high fc's 


RUN 1 
RUN 2 
RUN 3 


/o = 0.1 

fo = 0.1 
various 


none 

T>i = i9(k* -\k\)ip(k,t) 
P, = t^(V- 2 )V(x,0 


2\ = ^(V 2 )V(x,t) 
2\ = ^(V 2 )V(x,t) 
V h = iv h (V 2 )*il>(x,t) 



Table 1: Summary of the different numerical simulations performed. Here the hat operator 
: stands for the Fourier-space operator, Vh = 2 X 10~ 6 , Vi = 1 X 10~ 18 , /x = 1 x 10—4, 



fc* = 9Afc and the forcing is defined in (26) 



the numerical error in the algorithm is only due to the time-splitting 
The time step is always At = 0.5, chosen to be of the order of the shortest 
linear time 2ir /u max . 

Concerning the forcing and the dissipation, it appears convenient to con- 
trol T and T> directly in Fourier space to gain a wide inertial range. The 
forcing term acts to inject mass and energy in the system. In all simulations, 
J- modifies the first (linear) calculation half-step as follows, 



tp(k,t + At) 



h ma,x 

I 

"rnaxj 



^(k, t)e" l l k l 2Ai + Afoe 1 ^, k min < |k| < K 
^(k,t)e" l l k l 2Ai , |k| < fc min U |k| > K 

(26) 

where A ~ 1.62 x 10 3 . Thus the pumping add mass and energy in the 
ring region k min < |k| < /c max with the function <^(k, t) uniformly distributed 
in [0, 2tt) and statistically independent at each k-space and each time step. 
We choose the forcing at relatively large scales, with /c mm = 9Ak and fc max = 
10AA;. A dissipation at high wave numbers is added to halt the direct cascade 
and to prevent termalization effects. We find that an hyper-viscous term 
T>h = zz// l (V 2 ) 8, ?/'(x, t), where i//, = 2x 10~ 6 , is effective in absorbing the high-/c 
spectrum and in preventing the aliasing and the bottleneck effects. This term 
is also added to the linear operator to evaluate it efficiently in Fourier space. 
Finally, different types of dissipations at the large scales can be chosen. In 
this manuscript we report three different setups, briefly summarized in Table 
[TJ whose results are discussed in the following. 

4-1. RUN 1 - free condensate growth 

Here, we study the evolution of the GPE system, initially empty, without 
any dissipation at low k's. At every time step, the forcing term inputs mass 
AM(t) and energy AE(t), so that M and E start to grow. Even if the forc- 
ing is acting at a particular wave-number range, the nonlinear interactions 
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Figure 4: Wave-action spectrum ni£>(k,t) during initial stages of simulation. The time 
unit r refers to the linear time of the forcing scale. 



cause mass and energy transfers in Fourier space. We have chosen the forc- 
ing coefficient fo such that the transfers become efficient for the time step 
considered, namely it chosen to be sufficiently big to prevent sandpile effects 
characteristic to mesoscopic turbulence One can see in Fig. [4] that at 

early stages of the simulation the wave-action spectrum nmik, t) evolves and 
spreads over the wave-numbers space: the energy and the particles undergo 
a direct and an inverse cascades respectively. After an initial transient time, 
the linear energy density £n n (t) stops to grow, as clearly visible in Fig. |5} 
Linear energy weights at the high wave number part of the spectrum and this 
saturation is an evidence that its transfer to the small scales is now absorbed 
by the hyper-viscosity. On the contrary, as no dissipation is present at the 
large scales, the inverse particle cascade is not arrested and the nonlinear 
energy density £ n i(t) continues to increase. 

4.. 1.1. Spectra 

In Fig. [6] we present the ID wave-action spectrum at two different stages 
after the stabilization of the linear energy. The spectrum plotted with dashed 
line is taken at early stages, when the linear and the nonlinear energy densi- 
ties are comparable: at this point there is a good agreement with the WWT 
prediction, which is also plotted. The condensate component c (t) con- 
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Figure 5: Linear £u n (t) and nonlinear £ n i(t) energy densities during the RUN 1. 




Figure 6: Wave-action spectrum nio(k,t) at two different stages of RUN 1. The four- 
wave and the three- wave WWT predictions are also indicated by the lines with slopes — 1 
and —3/2 respectively. Inset: evolution of the condensate component co(t). 
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Figure 7: Dispersion relation at the final stage of RUN 1, when the condensate fraction 
is strong. The white dashed line is the Bogoliubov dispersion (only the upper branch). 
Inset: a zoom on the small k's zone to appreciate the shift due to the presence of the 
condensate. 



tinues to increase during the simulation, as shown in inset to Fig. [6j and a 
well defined series of peaks appears in the wave-action spectrum at the final 
stage (continuous line). Such behavior is a clear sign of three- wave interac- 
tions, as reported in [22! . Note that the late stage spectrum is consistent 
with the three-wave Zakharov-Sagdeev WWT prediction k~ 3 ^ 2 . The acous- 
tic three-wave regime is also confirmed by evaluating the dispersion relation 
by measuring the wavenumber- frequency Fourier transform. The result, pre- 
sented in Fig. [7J shows the presence of two branches (one for the Bogoliubov 
mode and another for its conjugate) shifted by the condensate velocity oscilla- 
tion Cq, in excellent agreement with the corresponding Bogoliubov dispersion 

The observed evolution of the spectrum can be summarized as follows, 
(a) At initial times the condensate wave amplitude ^k=o(£) is of the same 
order as the amplitude of other modes and the dynamics is well described by 
the four-wave WWT regime, (b) As the inverse cascade is not halted, mass 
accumulates at low k's and strong turbulence takes place in this transition, 
(c) Finally when the zero-mode becomes dominant over the fluctuations the 
wave turbulence is again weak and well described by the three- wave WWT. 



16 



~ 10 2 



t 10" 1 I 



10" 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 
Density p 



Figure 8: Probability density function of the density field p at different time: early stages 
(t a ) and final stage in RUN 1. A normal distribution with mean fi and standard 
deviation a are indicated. 



4-1.2. Condensation and density PDF 

To access further information about the statistics of the transition from 
the four-wave to the three-wave regime during the condensation process, in 
Fig. [8] we show the probability density function (PDF) of the density field p 
at initial and final stages. At early stages the density remains small with a lot 
of low density regions present, which indicates presence of numerous "ghost" 
(weakly nonlinear) vortices. In fact, in the ideal four- wave WWT dynamics, 
the wave field would be nearly Gaussian, and the respective density p would 
have an exponential (Rayleigh) PDF, which would be a straight line (with a 
negative slope) in Fig. [8j We see a significant deviation from the Rayleigh 
behavior, which means that even at an early stage (at time t a when the 
four-wave KZ spectrum is reported) the statistics already differed from the 
Gaussian. The development of a maximum on the PDF of p is a signature of 
the emerging condensate. As the inverse cascade is not halted the condensate 
density keeps growing in time. At final stages (at time tb when the three- 
wave KZ spectrum is reported) the density field shows a normal distribution 
behavior in the core of the PDF with a tail remaining at low density regions 
corresponding to vortices (which are now strongly nonlinear). Note that 
the mean value of p at this time (0.428) is much greater than the standard 
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Figure 9: Snapshots of iso-surfaces of low density region with threshold p t h r = 0.1 (the 
mean density is (p) ~ 0.43). The box is 1/2 x 1/2 x 1/4 the computational domain. The 
frames are taken approximately every 3r. 



deviation (0.055) which means that the condensate density is much stronger 
that the Boboliubov fluctuations about the mean density. This is a clear sign 
of the three- wave WWT. 



4-1.3. Vortices and velocity PDF 

Now let us focus on the vortex component represented by the low-p PDF 
tail at t = tb- By plotting the iso-surfaces with small density threshold 
{pthr = 0.1) only one quantum vortex is found in the computational box. We 
plot in Fig. [9] its evolution, showing only a part of the total box. The vortex 
has a ring shape and it propagates in the direction of the ring axis. The vortex 
core radius is consistent with the healing length estimate £ ~ 1.5 = 1.5Ax. 
Propagating Kelvin waves can be observed on the vortex line. We will argue 
below that these waves are crucial for understanding the 5/3-spectrum of the 
incompressible kinetic energy. 



In Fig. 10 we plot the PDF of the single velocity components (the data 



are normalized in order to compare different distributions). The velocity 
PDF appears to have a dominant Gaussian core, which is consistent with 
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Figure 10: Probability density function of the single velocity components at the final stage 
in RUN 1. Inset: zoom on the positive tail of the distribution in log-log coordinates. 



the three-wave WWT. In addition, the velocity PDF power-law tail with 
exponent —3 is another signature of the thin vortex lines, whose velocity 
field falls off inversely proportional to the distance from the line at short dis- 
tances. Such power-law velocity PDF's where observed in superfluid turbu- 
lence experimentally [IB] and numerically [29J and interpreted as an evidence 
of quantum vortices. 



If.. 1.4- Incompressible energy spectrum 

Spectra of the compressible and the incompressible kinetic energy at late 
time (t = tb) are presented in Fig. 11 The compressible part is dominant and 
has the same features as the ones already discussed for the wave-action spec- 
trum: it shows the peaks (harmonics of the forcing scale) and it follows the 
Zakharov-Sagdeev spectrum (with exponent —3/2 + 2 = 1/2). The incom- 
pressible spectrum formally coincides with the classical Kolmogorov 5/3-law. 
But we have seen in Fig. [9] that only one quantized vortex ring remains in 
the system at this time, so it is impossible for the Kolmogorov theory, devel- 
oped for the continuous classical vorticity fields, to be relevant in this case. 
Resolution to this paradox is suggested in [19]. In short, the —5/3 slope is 
produced by an energy cascade carried by Kelvin waves (seen in Fig. [9]). It 
turns out [27] that the Kelvin wave energy spectrum also has exponent —5/3 
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Figure 11: Compressible and incompressible kinetic energy spectra measured at final stage 
of simulation in RUN 1. A /j~ 5 / 3 showin classical Kolmogorov prediction is also plotted. 



and the coincidence with the Kolmogorov exponent is purely coincidental 
(the Kelvin- wave and the Kolmogorov spectra have different pre-factors) . 

4-2. RUN 2 - friction at large scales 

In the previous run no steady turbulent state has been reached because of 
the presence of the inverse cascade. To stay in the four-wave weak turbulence 
regime and to avoid the condensate growth, an effective friction term at large 
scales will now be added. This term is written directly in Fourier space and 
results in 

T>i = i»6{k* - |k|)$, (27) 

where 9 is the Heaviside step function, k* = 9Ak is the lowest forced wave- 
number and /i = 1 x 10~ 4 is a friction coefficient. This coefficient is optimally 
chosen to stop the inverse cascade without altering the direct energy cascade. 
With this choice both the nonlinear and the linear energy densities reach a 



constant value during the simulation, as plotted in Fig. 12 The final stage 



wave-action spectrum, presented in Fig. 13, agrees with the k~ l four- wave 
WWT prediction. The condensate growth, shown in the inset, is halted 
by the friction. Agreement with the four-wave WWT may seem surprising 



because, according to Fig. 12 , the nonlinear energy exceeds the linear one. 



Our explanation is that the nonlinear energy mostly resides in the condensate 
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Figure 12: Linear £n n (t) and nonlinear £ n i(t) energy density evolution in the presence of 
large scale friction (RUN 2). 



and forcing scales, whereas in the direct cascade range the modes are weakly 
nonlinear. 

We now turn our attention to the statistically steady state distributions 
in the physical space. The PDF of the density field p, not showed here, looks 
similar to early stage distribution presented in Fig. [SJ Here the average 
density is (p) ~ 8.31 • 10 -2 , which corresponds to healing length £ ~ 3.5Ax. 



The low density regions in the computational box are plotted in Fig. 14 A 



big fuzzy ring structure is present in the top of the figure. It is probably 
a single ring with large fluctuations, which create secondary small vortex 
loops near the main vortex. Besides this ring, other uniform bubble-like low 
density regions are present in the box. These are small scale ghost vortices 
which are weakly nonlinear and short-lived (their typical size is close to the 
resolution scale). Thus we see that even though the wave field is mostly 
random, the coherent vortex structures are also seen in this regime. 

The system in the fluid dynamics framework presents quite unexpected 
results, still remaining to be explained. The PDFs of the velocity compo- 
nents are plotted in Fig. 15 The PDF is isotropic and Gaussian in the core, 
but it has a power-law tail ~ fj~ 5 ' 3 . We still do not have a theoretical expli- 
cation for this exponent, although the power-law PDF tail with a different 
exponent (-1) was previously predicted for WWT in [30J. The compressible 
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Figure 13: Wave-action spectrum nm{k,t) at final stage of simulation in the presence of 
friction. The WWT four- wave prediction fc _1 is shown by the straight line. The inset 
shows the condensate evolution co(t). 




Figure 14: Low density region of the density field, threshold pthr — 0.015, in presence of 
friction (RUN 2). 
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Figure 15: Probability density function of velocity components in the presence of friction 
(RUN 2). Inset: zoom in log, log scale of the positive PDFs branch. 



and the incompressible kinetic energy spectra are presented in Fig. [T6j The 
compressible spectrum does not show peaks (as in RUN 1) indicating that 
this turbulence is not acoustic. The incompressible spectrum shows a power- 
law behavior with exponent close to 1/3. Again, no theoretical explanation 
could be proposed. 

4-3. RUN 3 - hypo-viscosity dissipation 

One can devise to stop the inverse cascade with a different type of dissi- 
pation at low fc's. For this, let us now choose a hypo-viscosity of the form 

Z>,=i(V- 2 )fy(x,t). (28) 

As this operator is singular in k = 0, we will separately remove, at each time 
step, the zero-mode in Fourier space. This choice still allows to reach steady 
turbulent states, but unexpectedly leads to different results with respect to 
RUN 2. In Fig. [17] we present various steady state spectra evaluated with 
different levels of the forcing / . What emerges clearly is that the spectra 
still have a power-law behavior in the inertial range, but it does not follow 
the WWT four- wave prediction k~ l . Instead, we see the k~ 2 critical balance 
prediction for the wide range of the forcing amplitudes. 
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Figure 16: Compressible and incompressible kinetic energy spectra in a steady state with 
friction (RUN 2). 




Figure 17: Wave-action spectra n\r>(k, t) at the final stage of the simulation in the presence 
of hypo- viscosity (RUN 3) for different forcing coefficients: fo = 0.05 (A), fo = 0.1 (B), 
fo = 0.5 (C), f = 1.0 (D), f = 2 (E), f = 3 (F). The fc" 1 WWT'and fc" 2 CB 
predictions are also shown. 
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Figure 18: Energy ratio ?/ = H n i/Hi in with respect to different forcing coefficients fo in 
the presence of hypo- viscosity (RUN 3). For informations about labels see caption in 

Fig. \n\ 



Our interpretation is the following. The hypo-viscosity causes an infrared 
bottleneck that alters the dynamics at the scales near the forcing: there the 
linear and the nonlinear energies become comparable. At this point the 
Fj0rtoft argument can no longer apply and the critical balance condition 
propagates into all the inertial range causing the observed k~ 2 spectra. This 
suggestion is corroborated by the measurement of the ratio between nonlinear 
and linear energies r\ = H n i/Hi in evaluated for various forcing coefficients fo 
illustrated in Fig. [18| From these results it is clear that, for a wide range of 
the forcing coefficients (almost two orders of magnitude), the energy ratio rj 
remains always order one. 

As in the previous runs, we look at the density field in the physical space 
to visualize vortices and other turbulent structures. For this we choose the 
case (D) when fo = 1.0 (all other cases look similar). Again, the PDF of 
the density is similar to the early stage PDF in Fig. [8j and this is natural 
because there is no condensate in the present system. The mean density is 
(p) ~ 4.31 ■ lO^ 1 and so £ ~ 1.5As. In Fig. [l9] we show the low density 
regions in the physical space with a threshold pthr = 0.05. This figure is 



qualitatively different from the previous ones (Fig. M and Fig. 14). Now 



very thin vortex structures fill completely the computational box and form 
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Figure 19: Low density regions with threshold of pthr — 0.05 are plotted at the final stage 
with the presence of hypo-viscosity with /o = 1.0 (RUN 3). 



a "vortex tangle" . 

The one point PDF's of the single velocity components are plotted in Fig. 



20 The PDF's, which show isotropy, are strongly non-Gaussian: the power- 
law behavior of v~ 3 is observed. This result is interesting because similar 
behavior was observed experimentally in jlH] and numerically in [29] and 
explained by presence of thin quantized vortex lines. We emphasize the v ~ 3 
PDF behavior is dominant and not present just in the tail as in the RUN 1. 
This is because we have much more strong vortex lines in the RUN 3 than 
in the RUN 1. Indeed, as the condensate fraction is removed by this type 
of dissipation, it is thus natural to think that the vortices fill the system at 
all scales and highly influence the velocity field. These vortex lines undergo 
frequent reconnections resulting in a sound emission. The incompressible 
kinetic energy spectra, for all the forcing amplitudes, are illustrated in Fig. 



21 No evidence of —5/3 law is found in this (critical balanced) regime. All 
spectra have a power-law behavior in the inertial range with exponent near 
— 1. A theoretical explanation of this exponent is still lacking. 
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Figure 20: PDF of single velocity components in the presence of hypo- viscosity with 
fo = 1.0 (RUN 3). Inset: zoom of the positive branch of the PDF in log-log coordinates. 
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Figure 21: Incompressible kinetic energy spectra for different forcing amplitude fo in 
presence of hypo- viscosity (RUN 3). For information about the labels see caption in Fig. 
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5. Conclusions 

In this paper, we have analyzed the turbulent states in the forced-dissipated 
3D GPE model by using the direct numerical simulations. Introduction of 
the forcing and the damping is aimed at achieving statistically stationary 
turbulent cascades. We have focused our attention on the direct energy cas- 
cade by introducing a pumping term at relatively large scales and by using 
hyper-viscosity at small scales. We have studied three regimes with different 
dampings at the scales larger than the forcing scale: no dissipation (RUN 
1), a friction (RUN 2) and a hypo- viscosity (RUN 3). 

RUN 1 is performed without any dissipation at large scales and so the 
inverse cascade, predicted by the WWT theory, causes condensation at the 
k = mode. This alters the four-wave dynamics and the system, after a 
strongly turbulent transient, becomes dominated by weak three- wave dynam- 
ics of acoustic fluctuations on background of a strong coherent condensate. 
The long-time evolution in this case is characterized by the appearance of 
a large quantum vortex ring in the numerical box. Kelvin waves propagate 
on this vortex ring, which causes the incompressible kinetic energy spectrum 
to follow the —5/3 law. We argue that the classical turbulence picture de- 
veloped for the continuous vorticity fields is inapplicable; the fact that the 
Kelvin wave turbulence has the same spectrum as the classical Kolmogorov 
spectrum is coincidental, as it arises from completely different physical pro- 
cesses j2?l S9]. The vortex ring in this regime coexists with the random 
acoustic waves engaged in the three-wave interactions. The vortex shows 
up as a v ~ 3 tail on the velocity PDF while the random waves make up the 
Gaussian core of this PDF. 

In RUN 2 and RUN 3 we have introduced a dissipation term at large 
scales in order to stop the inverse cascade and to reach statistically steady 
states. The characteristics of these states depend strongly on the choice 
of the low-fc damping. If a friction is introduced (RUN 2), the growth 
of the condensate is halted and the wave-action spectrum follows the four- 
wave WWT prediction. A large fuzzy vortex ring surrounded by small ghost 
vortices appears in the final stage of the computation. If the dissipation is 
an hypo- viscosity (RUN 3), the final steady spectra, evaluated for a wide 
range of forcing coefficient, agree with the critical balance prediction. In 
this regime the computational box appears to be filled by a vortex tangle - 
a chaotic set of strongly nonlinear vortex lines. The velocity PDF exhibits, 
both in the core and on the tails, a power-law behavior v~ 3 characteristic to 
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such vortex lines. 

In Summary, our numerical results clearly show that the turbulent state 
in the direct cascade range is strongly affected by the choice of damping in 
the inverse cascade range. Most realistic configuration for the existing BEC 
experiments is configuration of RUN 1. Indeed, while the hyper- viscosity 
can be physically understood as an evaporative cooling mechanism, no large- 
scale damping mechanisms have ever been proposed. On the other hand, to 
study the nontrivial QT states predicted by the RUN 2 and the RUN 3 
of the present paper, it would be interesting to explore possibilities to damp 
the lowest-momentum modes in BEC experiments. 
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Appendix A. Signs of the fluxes 

The KZ solutions carry constant fluxes of the conserved quantities over 
the turbulent scales. The GPE model has two conserved quantities: the mass 
(particles) and the energy. Here, we will find the direction of the fluxes cor- 
responding to the KZ solutions in a very simple way, avoiding computing the 
flux based on a complicated relation resulting from the the kinetic equation, 
as it was done in [50] and [33J. 

Consider Fig. A.22| where we qualitatively plot the flux on a generic 



power-law spectrum n{k) = ck~ a (which is note necessarily a steady solu- 
tion) as a function of the exponent a. For very sharp spectra, a>l, both 
fluxes must be positive. To see that one can think of a narrow-band spectrum: 
it must spread which corresponds to the positive fluxes on the negative slope 
side. As the fluxes are continuous functions of a, we can determine their di- 
rections based on the zero-crossing points, i.e. the RJ and the KZ exponents. 
Of course, at the two thermodynamic RJ solutions, corresponding to a = 
and a = 1, both fluxes must be zero. At the KZ solutions, only one of the 
fluxes turns into zero. Namely, on the energy cascade KZ, the mass flux is 
null and vice versa. Thus, the functions rj(a) and e(a) qualitatively behave 



29 







a 



Figure A. 22: Qualitative direction behaviors of energy and particles fluxes with respect to 
the power-law exponent a of KZ solution n(k) = ck a . 



as shown in Fig. A. 22 It is then clear that the energy undergoes a direct 



cascade while the mass cascade inversely. 

Appendix B. Locality of interactions 

Lets test the locality of interactions in the constant-flux states. This imply 



checkig that the collision integral in the kinetic equation (13) converges on 
the KZ spectra n(k) = ck~ a . As these solutions are scale-invariant, the 
integral is easily written in the u space as 

X5(ui +U 2 -U 3 - UJ 4 )f{Ul, U 2 , U 3 , U 4 )du234, (B.l) 

where term 

f( , min yu[, yfo^ y/uj, y/UT^ 

f(u 1 ,u 2 ,uJ 3 ,uj i ) = -j (B.2) 

takes into account the 3D average of 5(ki + k 2 — k 3 — k 4 ) over the solid angles 
Qi, Q, 2 , fl 3 and ^4 , see [3j [33l EH E2] for details. In this coordinate system 
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the frequency 8- function can be easily used, for example as lo 2 = w 3 + cj 4 — U\ . 
The integral presents singularities for integration over at zero and infinity 
Note that is the same is true for U4 integration as the integral is symmetric 
with respect to (3, 4) — > (4, 3). By Taylor expanding the integrand up to the 
leading order, one gets 



_ a + 2 

lim ~ / (X> 3 2 du3 and 



W3— »oo 



lim ~ / uj 3 2 du^. 

us— > 



(B.3) 



The locality holds for < a < 3 which is true for the inverse particles cascade 
(a = 7/3) but not for the direct energy cascade (a = 3). Nevertheless the 
divergence in the latter case is marginal and the solution can be corrected 
by a logarithmic factor |33j. 

Appendix C. The Fj0rtoft argument 

We will present here a new version of the argument Fj0rtoft argument 
which is formulated for the conservative systems (no forcing or dissipation). 
Let us introduce the mass and the energy centroids as 

Jkn(k)dk 

= M ^ ' 

jk£{k)dk fk 3 n(k)dk 
K E =~ ]^ = ~ • ( C - lb ) 

In the latter expression we have assumed (this is essential for the Fj0rtoft 
argument) that the nonlinear energy is negligible with respect to linear one. 
In the following will use the Cauchy-Schwartz inequality 

V 2 / r \ l l 2 



f(k)g(k)dk < (J f 2 (k)dkj (J g\k)dk ) (C.2) 
By splitting the energy integrand in two parts we get 

E = [ k 2 n(k)dk = [ (k l/2 n 1/2 ) x (k 3/2 n 1/2 )dk 



1/2 x r V 1/2 

< I / kn(k)dk ) / k 3 n(k)dk = y/K M MK E E (C.3) 
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and so we have 

K U K E > ^. (C.4) 

This inequality means that if the particle centroid moves to low wave numbers 
(inverse cascade) the energy centroid must move to high wave numbers (direct 
cascade). 

We now evaluate 

K M M = [ kn{k)dk = [ (kn 1/2 ) x {n 1/2 )dk 



1/2 , f \ 1/2 

< ( / k 2 n(k)dk ) / n(k)dk = VME (C.5) 



which gives 



K M < \[EjM. (C.6) 



This inequality means that the mass centroid can either stay where it is 
initially, or move to the large scales (inverse cascade), but it cannot cascade 
to the small scales. 



Combining the inequalities (C.4) and (C.6), we get 



K E > s/EjM. (C.7) 

This inequality means that if the energy centroid can either stay where it 
was initially, or move to the smaller scales (direct cascade), but it cannot 
cascade to the larger scales. 
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